Universal Conductance Fluctuations in Mesoscopic Systems with Superconducting 

Leads: Beyond the Andreev Approximation 
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We report our investigation of the sample to sample fluctuation in transport properties of phase co- 
herent normal metal-superconductor hybrid systems. Extensive numerical simulations were carried 
out for quasi-one dimensional and two dimensional systems in both square lattice (Fermi electron) as 
well as honeycomb lattice (Dirac electron). Our results show that when the Fermi energy is within 
the superconducting energy gap A, the Andreev conductance fluctuation exhibits a universal value 
(UCF) which is approximately two times larger than that in the normal systems. According to the 
random matrix theory, the electron-hole degeneracy (ehD) in the Andreev reflections (AR) plays an 
important role in classifying UCF. Our results confirm this. We found that in the diffusive regime 
there are two UCF plateaus, one corresponds to the complete electron-hole symmetry (with ehD) 
class and the other to conventional electron- hole conversion (ehD broken). In addition, we have 
studied the Andreev conductance distribution and found that for the fixed average conductance (G) 
the Andreev conductance distribution is a universal function that depends only on the ehD. In the 
localized regime, our results show that ehD continues to serve as an indicator for different universal 
classes. Finally, if normal transport is present, i.e., Fermi energy is beyond energy gap A, the AR 
is suppressed drastically in the localized regime by the disorder and the ehD becomes irrelevant. As 
a result, the conductance distribution is that same as that of normal systems. 

PACS numbers: 72.80.Vp, 74.45. +c, 73.23.-b, 68.65.Pq 
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I. INTRODUCTION 

It is well known that quantum interference leads to 
significant sample-to-sample fluctuations in the conduc- 
tance at low temperatures. These fluctuations can be 
observed in a single sample as a function of external pa- 
rameters such as the magnetic field since the variation of 
magnetic field has a similar effect on the interference pat- 
tern as the variation in impurity configuration. One of 
the fundamental problems of mesoscopic physics is to un- 
derstand the statistical distribution of the conductance 
in disordered systems 2 . - — . 

It has been established that in the diffusive regime, 
the conductance of any metallic sample fluctuates as a 
function of chemical potential, impurity configuration 
(or magnetic field) with a universal conductance fluc- 
tuation (UCF) that depends only on the dimensional- 
ity and the symmetry of the system^ The UCF is given 
by Var(G/G ) = 2/(16/3), 2/(15/3), 3/(16/3), 5/(17/3) for 
quantum dot (QD), quasi-one dimension (ID), two di- 
mensions (2D) square and three dimensions cubic sample 
with Go = 2e 2 /h. Here the index /3 corresponds to circu- 
lar orthogonal ensemble (COE) when the time-reversal 
and spin-rotation symmetries are present (J3 — 1), circu- 
lar unitary ensemble (CUE) if time-reversal symmetry is 
broken (/3 = 2) and circular symplectic ensemble (CSE) if 
the spin-rotation symmetry is broken while time-reversal 
symmetry is maintained (/3 = 4), respectively £ In the 
crossover regime from diffusive to localized regimes, the 
conductance distribution was found to be a universal 
function that depends only on the average conductance 
for quasi-lD, 2D, and QD mesoscopic systems and for 



(3 = 1,2,4£I In the localized regime, the conductance 
distribution seems to be independent of dimensionality 
and ensemble symmetry^ 

In the presence of a superconducting lead, using ran- 
dom matrix theory (RMT), the conductance fluctuations 
in the mesoscopic normal and superconductor hybrid sys- 
tems have been studied in the diffusive regime for quasi- 
lD systems^ - — and QD system^ It was found that the 
UCF in a COE system shows approximately a twofold 
increase over the normal systems, i.e., rms(GArs) — 
2 rms(GAr)ji&ii Different from the normal conductor, in 
the presence of superconducting lead, electron- hole de- 
generacy (ehD) plays a similar role of "symmetry" . UCF 
assumes different value depending on whether ehD is 
broken or not. According to RMT, 3 the Andreev con- 
ductance fluctuation rms(GAT5) = \/4-3 rms(GAr) (with 
ehD) and rms(G_/vs) = rms(G.iv) (with ehD broken) 
were predicted. Up to now, however, most of the inves- 
tigations on Andreev conductance fluctuation have been 
done for systems with ehD (e = 0) and low energy regime 
(A -c E c where E c is Thouless energy). There is not yet 
a numerical study on the NS hybrid system where ehD is 
broken. In fact, for the existing studies on the NS hybrid 
system with ehD, there is no consensus on the theoreti- 
cal predicted value of rms(GArs). Specifically, concerning 
the increase factor ag in "rms(Gjvs) == a o rms(GAr)", a 
diagrammatic theory predicted -\/6j^ and a numerical 
calculation using tight binding model gave «o = V4, 10 
and the random matrix theory indicated cvo = V4.3^ and 

Recently, graphene based normal-metal- 
superconductor (GNS) systems were intensively 
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studied because good contacts between the super- 
conductor electrodes and graphene have been realized 
experimentally. 13,14 In a conventional (quadratic energy 
dispersion relation) normal-metal-superconductor (CNS) 
system, the usual Andreev reflection (AR) occurs^ For 
GNS systems, AR can be either intravalley or inter- 
valley, which are called Andreev retroreflection (ARR) 
and specular Andreev reflection (SAR), respectively^ 
When the excitation energy 2e is smaller than that of 
incident energy relative to Dirac Point Ep — Eq, ARR 
happens, otherwise SAR occurs. At the transition point 
2e = Ep — Eq between ARR and SAR, the reflection 
angle 9 (measured relative to the NS junction normal) 
jumps from +90° to —90°^ the shot noise vanishes and 
the Fano factor has a universal value.— In general, SAR 
differs from ARR or conventional AR (CAR) where an 
extra phase n which can be observed in the quantum 
interference of the two SAR reflections.— 

So far most of investigations on UCF focus on the 
Fermi electrons (quadratic dispersion relation) with the 
zero or low energy and less attention is paid on the Dirac 
electrons. In addition there is no numerical work report- 
ing Andreev conductance fluctuation when ehD is bro- 
ken. It would be interesting to ask the following ques- 
tions. What happens to UCF for GNS systems? Is it the 
same as that in CNS systems? Is there any difference 
between ARR and SAR? Which theoretically predicted 
value of UCF for the quasi-lD CNS system (with ehD) 
is favored? What happened when ehD is broken? What 
about the conductance distributions in these systems? It 
is the purpose of this paper to address these questions. 

In this paper, using the tight-binding model, we carry 
out a theoretical study on the sample to sample fluc- 
tuation in transport properties of phase coherent sys- 
tems with normal metal-superconductor heterojunction. 
In view of the possible difference among CAR, ARR 
and SAR, we consider both the CNS systems using the 
square lattice and GNS system using the honeycomb lat- 
tice. Extensive numerical simulations on quasi- ID and 
2D systems in the presence of a superconducting lead 
show that when the Fermi energy is within the super- 
conducting gap Ep < A, UCF roughly doubles the value 
in the absence of the superconducting lead. This is the 
case for both CAR in CNS system and ARR and SAR 
in GNS system. So there is no distinct difference be- 
tween ARR and SAR. Besides, concerning ehD in the 
NS hybrid system, new universal classes are present in 
agreement with the prediction of RMT4 Two plateaus 
of UCF were found in our numerical results, one corre- 
sponds to the complete electron-hole symmetry^ class 
(with ehD) and the other to conventional electron-hole 
conversion (with ehD broken). It was found that the 
case of "ehD broken" decreases the value of UCF, again 
in agreement with the theoretical analysis^ Specifically, 
in the quasi-lD systems, rms(Givs)/rms(Gjv) for both 
Fermi and Dirac electrons is 2.07 ± 0.04 that is close to 
V 4.3 when ehD is present while when ehD is broken it is 
1.99±0.08 that is close to \fi. For 2D systems, when ehD 
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FIG. 1: (Color online) Schetch of CNS [panel (a)] and GNS 
[panel (b)] system, in which ideal superconducting lead (left, 
orange), normal lead (right, blue) and disordered normal scat- 
tering region (shadowed blue region) are concluded. 



is present, rms(Gjv,s)/rms(GAr) is 1.91 ± 0.07 for Fermi 
electrons and 1.96 ± 0.07 for Dirac electrons while it is 
1.82 ± 0.08 when ehD is broken for both Fermi electrons 
and Dirac electrons. Furthermore, the different conduc- 
tance distributions -P(G) for the fixed average conduc- 
tance (G) also indicate this new symmetry class in local- 
ized regime. We also point out that the new universality 
class due to the ehD is quite different from the conven- 
tional ensemble symmetries. It was shown numerically 
that the conductance distribution P(G) in the deep lo- 
calized regime for normal systems is a universal function 
which depends only on the average conductance (G) but 
not on the Fermi energies as well as other parameters. 7 
In addition, it does not seem to depend on the ensemble 
symmetry and dimensionality of the system. In the pres- 
ence of the superconducting lead, our numerical results 
for 2D systems with (3 = 1 show that the conductance 
distribution is still an universal function that depends 
only the average conductance (G). Different from nor- 
mal system, however, it depends on whether the system 
has the ehD. Finally, when Ep is above A, normal trans- 
port is present. We found that the AR is suppressed 
by the disorder especially in the localized regime where 
normal transmission dominates transport processes. In 
this case, the ehD is irrelevant and the same universal 
conductance distribution is found as that in the normal 
systems in the localized regime. 

The rest of the paper is organized as follows. In Sec. II, 
with the tight-binding representation, the model system 
including central disordered region and attached ideal 
normal lead and superconducting lead is introduced. The 
formalisms for calculating the conductance and fluctua- 
tion of conductance arc then derived. Sec. Ill gives nu- 
merical results along with detailed discussions. Finally, 
a brief summary is presented in Sec. IV. 



II. MODEL AND HAMILTONIAN 

The scattering theory of electronic conduction is de- 
veloped by Landauer^i Imry^ and Biittiker<22 It pro- 
vides a complete description of quantum transport in the 
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system without electron-electron interactions. A meso- 
scopic conductor can be modeled by a phase-coherent 
disordered region connected by ideal leads (without dis- 
order) to two electron reservoirs (normal metal or super- 
conductor) , which are in equilibrium at zero temperature 
with fixed electrochemical potential (or Fermi energy) 
E F . Here we assume that the central scattering region is 
normal region, the same as the right normal lead. Then 
the total system Hamiltonian 



H = H S + H 



N 



(1) 



where H$, Hn and H F are the Hamiltonian of supercon- 
ducting lead (orange region in Fig.l), semi-infinite nor- 
mal ribbon (blue region in Fig.l) and tunneling between 
the normal region and superconducting terminal, respec- 
tively. 

Two kinds of structure were considered in this paper: 
the structure with quadratic energy dispersion [square 
lattice, Fig. 1(a)] and structure with conical energy spec- 
trum [honeycomb lattice, Fig. 1(b)]. In the absence 
of the superconductor, the whole system Hq including 
Hs, Hn and He can be written in the tight-binding 
representation ; 24 ' 25 
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where i = (i x ,i y ) is the index of the discrete square lat- 
tice or honeycomb lattice site which is arranged as in 
inset of Fig.l. Here a\ and a\ are the annihilation and 
creation operators at the discrete site i. Eq is the con- 
stant on-site energy. In the square lattice, Eg is center of 
energy band, and in honeycomb lattice, Eq is the energy 
reference point (the Dirac Point). Sei is random on-site 
potential which is nonzero only in the center region to 
simulate the disordered scattering region. Here Sei is 
uniformly distributed with Sei = [— w/2, w/2] where w is 
disorder strength. The data for fluctuations are obtained 
by averaging over up to 10,000 disorder configurations 
and the data for distribution are obtained over 1,000,000 
disorder configurations. The second term in Eq.(2) is 
the nearest neighbor hopping with hopping elements t 
and "<>" denotes the sum over the nearest sites. 

Due to the superconductor, it is convenient to write the 
Hamiltonian H in the Nambu representation^. In this 
representation the Fermi energy of the right normal lead 
in equilibrium (at zero bias) is set to be the supercon- 
ductor condensate. It is conventionally set to zero. As 
a result, the spin up electrons and the spin down holes 
have the positive and negative energy, respectively. Tak- 
ing this into account the Hamiltonian ^ is cross multi- 
plied by spin representation. Hn and Ht in Eq.JT]) can 
be rewritten as H n/t ® o z and 



H s = 



H Q .s A 
A* -H , s 



(3) 



where A = Ae lip is the energy gap or the pair potential of 
the semi-infinite superconducting lead. Here we can as- 



sume A = A to be a real parameter by selecting a special 
phase of the superconductor lead in our calculation^ 

In the calculation, for simplicity we set external voltage 
in the normal and superconducting terminal as Vn = V , 
Vs = 0. The current flowing from the normal lead can 
be calculated from the Landauer-Biittiker formula^ 



Jn — J 

T e/h 
N 



N 



±- 



rh 



dE 
2^ 



{T e/h (E)[f ± (E) - f Q (E)} 
T A {E)[.f ± {E) - f T (E)]} (4) 



where e is the electron charge, fo{E) — [ e B / fcisr + l] 
is the Fermi distribution in the superconducting lead, 
f±(E) = [ e ( £= FeViv)/fc B r + are the Fermi distribu- 
tion functions in the normal terminal for the electrons 
and holes, respectively. T e / h is the transmission coef- 
ficient that the particles incident from superconducting 
lead traverse to the normal terminal as electrons/holes 
and Ta is AR coefficient representing the reflection prob- 
ability that the incident electrons from the normal ter- 
minal are reflected as holes or vice versa. Note that the 
two processes are symmetric and have the same AR co- 
efficient Ta- T e / h and Ta are calculated from 

T e = Tr{rf t [G"T s G a ] tt }, T h = Tr^CT 5 G a ] u }} 
Ta = Tr[r^.G* , rf, G%] = Tr[rjX GT*r^.G*, ] (5) 



the line- width function T N / S {E) = i[E r N/s {E) - 
E r l s (E)]. The Green's function G r {E) = [G a {E)}^ = 

[EI - H c - T, r N (E) - ^(E)}- 1 where H c is Hamil- 
tonian matrix of the central scattering region and I is 
the unit matrix with the same dimension as that of 
He, S[ =Ar g is the matrix of retarded self-energy from 
the normal/superconducting lead with the only nonzero 
elements in the sub-block that are neighbor of normal 
or superconducting lead. The self-energy is calculat- 
ing according to £[ = H C ig r Hi C where H C i [Hie) is 
the coupling from central region (leads) to leads (cen- 
tral region) and g r is the surface retarded Green's func- 
tion of semi-infinite lead which can be calculated using a 
transfer matrix method* 2 ^ Due to electron-hole symme- 
try, T e (E) = T h (-E) and T A {E) = T A {-E), which leads 
to Jn — 2Jn — — 2</jy-. 

At zero temperature limit, the energy dependent con- 
ductance can be expressed as: 

G N s(E F ) = d(J%-J N )/dV 



— {[T e (E F ) + T h {~E F )] + 4T A (E F )} 



2e 2 



[T e/h (±E F ) + 2T A (E F )] 



(6) 



When the incident energy E F < A, there is no nor- 
mal quasi-particle transport T e j h — and only AR con- 
tributes to conductance G. We will focus mainly on this 
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quantity in this paper. In this case, the conductance 
fluctuation defined as rms(G) = ([G — (G)] 2 ) becomes 

rms(G) = —j^-yj (T%) ~ (Ta) 2 (7) 

where (...) denotes averaging over an ensemble of sam- 
ples with different disorder configurations of the same 
strength w. When Ep is beyond superconducting energy 
gap A, normal transmission T e /^ is present, conductance 
variance now consists of three components: (1). the An- 
dreev related fluctuation Var(G)Andr from AR coefficient 
Ta- (2). the normal fluctuation Var(G)Norm from nor- 
mal transmission coefficient T e / h . (3). the cross term 
Var(G) cross . They are expressed as 

Var(G) = Var(G) A ndr + Var(G) Norm + Var(G) cross 
= [j? [((4«a) 2 ) + <(«V) 2 ) + 8 (ST A 5T N }} (8) 
where 6T A =T A - (T A ), 5T N = (T e + T h ) - (T e + T h ). 

III. RESULTS AND DISCUSSION 

In the numerical calculations, the energy is measured 
in the unit of the nearest coupling elements t. For the 
square lattice, t — 2t ", ^ with m* the effective electron 
mass and a the lattice constant. For the honeycomb lat- 
tice, t = -^fovp with the carbon-carbon distance b = 
0.142nm and the Fermi velocity vf = 0.89 x 10 6 ms -1 . 
The size of the scattering region N x M is described by 
integer N and M corresponding to the width and length, 
respectively. For example in Fig.l, the width W = Na 
with N = 3, the length L = Ma with M — 5 in the panel 
(a), and the width W = N x 3b with N = 3, the length 
L = M x \/36 with M = 7 in the panel (b). 

As documented in the literature, in order to get the 
saturated UCF plateaus j£i>2ii the number of transmission 
channels for incoming electron should be large enough in 
the numerical calculation. We denote N c as the chain 
number which determines directly the number of chan- 
nels. N c is defined in the following way: take Fig.l as an 
example, in panel (a), N c = 3 and N c = 6 in panel (b). In 
2D systems we set N c = 40 and 60. For quasi ID systems 
we use only N c — 40 because it is more computational 
demanding than 2D systems. To get a larger channel 
number, the incident energy Ep should be set away from 
the bottom of energy band Ef,. In the square lattice, 
to mimic the parabolic energy spectrum for Fermi elec- 
trons, the constraints for incident energy Ep < Ef, + 2t 
and Andreev reflected energy —Ep < Ef, + 2t are needed. 
While for Dirac electrons in the honeycomb lattice, the 
absolute value of relative incident energy (to the Dirac 
point Eq) \Ep — Eq\ < t and relative Andreev reflected 
energy | — Ep — E \ <t are set. 

In TableJH we list all the parameters used in the fol- 
lowing calculations including the incident energy Ep, the 
superconducting gap A and the on-site energy Eq which 



TABLE I: The parameter Ef, Eq, A used in the square lat- 
tice model and honeycomb model. The different columns are 
corresponding to the different transport processes denoted by 
'AR', 'ARR', 'SAR' ,'NT and so on. Here, 'AR' is for the 
pure conventional AR assisted tunneling processes (only con- 
ventional AR exists) in square lattice. 'ARR' and 'SAR' de- 
notes the pure ARR assisted process and pure SAR assisted 
process in honeycomb lattice, respectively. 'NT' is for the 
transport beyond the superconducting Gap where NT can 
also contribute to the transport processes. 



sq 



AR E F 


E 


A 


AR 


Ef 


E 


A 


NT 


E F 


E 


A 


1 


2.1 


0.1 


4 


0.2 


2.2 


0.3 


1 


0.2 


2.2 


0.3 


2 


2.3 


0.1 


5 


0.3 


2.3 


0.4 


2 


0.3 


2.3 


0.4 


3 


2.4 


0.1 


6 


0.4 


2.4 


0.5 


3 


0.4 


2.4 


0.5 


he 


ARR E F 


E 


A 


SAR E F 


E 


A 


NT 


E F 


E 


A 


1 


0.6 


0.1 


1 


0.6 


0.0 


0.7 


1 


0.7 





0.5 


2 


0.7 


0.1 


2 


0.6 


0.1 


0.7 


2 


0.7 


0.1 


0.5 


3 


0.8 


0.1 


3 


0.7 


0.0 


0.8 


3 


0.7 





0.3 


4 


0.9 


0.1 


4 


0.7 


0.1 


0.8 


4 


0.7 


0.1 


0.3 


5 0.1 


0.7 


0.2 










5 


0.7 





0.1 


6 0.1 


0.8 


0.2 










6 


0.7 


0.1 


0.1 



is the center of energy band for the square lattice and 
the Dirac Point for the honeycomb lattice. From these 
parameters in the clean system with NS heterojunction, 
we can easily calculate the channel number for electron 
or hole, AR coefficient Ta and the normal transmission 
coefficient of electron or hole T se / sh for Fermi energy be- 
yond superconducting gap. At the same time, we can also 
get the normal transmission coefficient T ne /„/j in normal 
system without the superconducting lead. 



A. Conductance fluctuation and conductance 
distribution in the diffusive regime 

We first examine conductance fluctuations in the dif- 
fusive regime. In our calculation the size of 2D square 
lattice is set to be 40 x 40 for N c = 40 and 60 x 60 for 
N c = 60. The size of 2D honeycomb lattice is chosen to 
be 20 x 35 for N c = 40 and 30 x 52 for N c = 60. For quasi- 
1D systems, the size is chosen to be 40 x 1000 in square 
lattice and 20 x 500 in honeycomb lattice with N c = 40. 
In Fig[2]and[3l [4]and[5j we plot conductance fluctuations 
rms(G) vs the average conductance (G) in 2D square lat- 
tice, quasi-lD square lattice, 2D honeycomb lattice and 
quasi-lD honeycomb lattice, respectively. Each point in 
the figure is obtained by averaging over 10,000 configu- 
rations. Different parameters used in all figures are tab- 
ulated in Table HI 

From FigJ5][Sl we see following general behaviors. (1). 
in the localized regime where (G) < 1, all the curves 
collapse into a single curve indicating the universal be- 
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FIG. 2: (Color online) panel(a) and panel(b): in the presence 
of superconducting lead, conductance fluctuation rms(G) vs 
average conductance (G) in the square lattice for N c = 40 and 
N c = 60, respectively. The symbol is labeled in first column 
in Tab|T]for the square lattice case denoted by 'sq'. The red 
dotted lines indicate two plateaus in the values (with the unit 
of 0.866e 2 /h) of 1.94 ± 0.03 and 1.85 ± 0.05 in panel (a) and 
1.98 ± 0.02 and 1.89 ± 0.03 in panel (b). For comparison, 
corresponding to panel (a) and panel (b), in panel (c) and 
(d), we plot rms(G) vs (G) in the absence of superconducting 
lead, i.e., A = 0, respectively. The plateaus in the values 
of 1.02 ± 0.02 and 1.04 ± 0.02 in the unit of 0.866e 2 //i are 
indicated in panel (c) and panel (d). The system size: N c = 
40 corresponds to width W = 40a, considering the square 
shape sample, we set L = 40a. For N c — 60, we have W = 
60a, L = 60a. 



FIG. 4: (Color online) rms(G) contributed by ARR [panel(a)] 
and SAR [panel(b)] vs (G) in 2D honeycomb lattice for 
N c = 40 [open symbols] and iV c = 60 [symbols with '-']. The 
symbols are labeled in second column in Tab|T]for the honey- 
comb lattice case denoted by 'he'. The red dotted lines indi- 
cate two plateaus with the values of 1.96±0.03 and 1.82±0.02 
in the unit of 0.866e 2 //i in panel (a) and a single plateau with 
the value of 1.82±0.02 in panel (b). Panel (c) and (d): rms(G) 
contributed by normal quasi-particle transmission T ne vs (G) 
in case of A = 0, corresponding to panel (a) and (b), respec- 
tively. The plateaus in the values of 1.00 ± 0.02 in the unit 
of 0.866e 2 //i are indicated in panel (c) and panel (d). The 
system size: for N c = 40 its width is equal to W = 606; 
considering the square shape sample, the length L — 35\/36. 
Similarly, for N c = 40 its width is W = 906, L = 52^6. 




10 20 30 40 



<G> (e /h) 

FIG. 3: (Color online) Same as Fig[2] except the model is 
quasi-lD square lattice with chain number N c = 40. The 
system size: width W = 40a, length L = 1000a. In the unit 
of 0.73e 2 //i, two plateaus with the values of 2.01 ± 0.02 and 
1.93 ± 0.03 in panel (a) and a single plateau in the value of 
0.97 ± 0.01 is indicated in panel (b). 
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FIG. 5: (Color online) Same as Fig[2] except the model is 
quasi-lD honeycomb lattice with chain number N c = 40. The 
system size: width W = 606, length L = 500^6. The red 
dotted lines indicate two plateaus with the values of 2.03±0.03 
and 1.93 ± 0.03 in the unit of 0.73e 2 /h in panel (a), a single 
plateau with the value of 1.94 ± 0.04 in panel (b), the value 
of 0.98 ± 0.02 in panel (c) and panel (d). 
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havior of the conductance distribution function^ (2). in 
the diffusive regime where (G) > 1, there is a plateau 
region for rms(G) where the fluctuation is nearly inde- 
pendent of average conductance (G) and other system 
parameters. This is the regime for the universal conduc- 
tance fluctuation. The plateau value [labeled by red dot- 
ted line in top panels] is approximately twice the value 
of the known UCF values rms(G) = b.8,66e 2 /h for 2D 
system and Q.73e 2 /h for quasi-lD system [labeled by red 
dotted line in bottom panels]. This doubling seems to 
be true for both Fermi electrons (square lattice) [Figf2] 
and FigJS] and Dirac electrons (graphene system) [FigfJ] 
and Figj5]. (3). There are two separate UCF plateaus 
for the AR assisted transport processes in the CNS sys- 
tem [FigfSJa), Fig[3ta)] and the ARR assisted transport 
processes in GNS system [FigJ2fa),FigJ5fa)], while for the 
SAR assisted transport processes [panel (b)] there is only 
one UCF plateau. It appears that this difference in UCF 
can be used to distinguish ARR and SAR. However, it 
turns out to be incorrect when considering the ehD "sym- 
metry". In FigfJIa) and FigJSJa), Andreev conductance 
fluctuations corresponding to Ep = (with ehD) and 
Ep 7^ (ehD broken) from diffusive regime all the way 
to localized regime are plotted and two UCF plateaus 
associated to ehD "symmetry" are then indicated. For 
SAR in graphene systems, we have \Ep\ > \Eq\ ^ (ehD 
broken). FigfJIb) and FigfSJb) then show only one UCF 
plateau. (4). Denoting the increase factor a through 
the relation rms(GArs) = ao rms(GAr) [Gat is shown in 
bottom panels in Figf2][5] in the plateau region in dif- 
fusive regime, it is very different for square 2D system 
and quasi- ID system and slightly different for Fermion 
electrons and Dirac electrons. Specifically, our results 
for the quasi- ID systems for both Fermi electrons and 
Dirac electron is as follows: (a), when ehD is present 
rms(Gjvs)/rms(GAr) is 2.07 ± 0.04 that is very close to 
\/4T3. (b). when ehD is broken it is 1.99 ± 0.08 that 
is close to \pi. For 2D systems, when ehD is present, 
rms(GAT5)/rms(GAr) is 1.91 ± 0.07 for Fermi electrons 
and 1.96 ± 0.07 for Dirac electrons. When ehD is bro- 
ken it is 1.82 ± 0.08 for both Fermi electrons and Dirac 
electrons. (5). For larger (G) (in the ballistic regime) 
the conductance fluctuation falls down quickly to zero. 
This is because the number of conducting channels N c is 
nnite£i2£. The width of plateau region is longer with a 
larger N c . In the limit of the infinite N c , the plateaus of 
conductance fluctuation will extend to infinite. 

We now take a closer look at each figure discussed 
above. In the top panels of Fig we can see that 
all curves of rms(G) vs (G) collapse into universal curves 
that are slightly separated in the region of 1 < (G) < 10. 
To make the discussion of separate UCF plateaus quan- 
titative, we plot rms(G) vs small (G) ((G) < 10) in 
FigUJa), (b), (c) and (d) corresponding to FigO FigEl 
Figf4] and FigJS] In Fig|6l we clearly see two separate 
UCF in the regime where 1 < (G) < 10. For each UCF 
plateau, the conductance fluctuation rms(G) vs average 
conductance (G) is a universal function, i.e., it is inde- 
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FIG. 6: (Color online) Corresponding to Fig[2] FigEJ Figg] 
and FigO rms(G) vs small (G) (< 10) are plotted in panel(a), 
(b), (c) and (d), respectively. 
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FIG. 7: (Color online) In the diffusive regime, corresponding 
to eight selected parameters with \Ef\ < A from Tab|H the 
conductance distribution P(G) obtained from 1,000,000 con- 
figurations is plotted for the fixed (G) ~ 3 in square [panel 
(a)] and honeycomb lattices [panel (b)]. 



pendent of system parameters such as Ep, Eq, A, system 
size and so on and depends only on (G). In fact, not 
only the rms(G) (the second moment), the third, forth, 
and higher moments are universal function of (G). 
This means that the conductance distribution P(G) is 
a universal function that depends only on the average 
conductance (G) in addition to the symmetry and di- 
mensionality of the system. 

To demonstrate the conductance distribution has two 
different universalities, we plot in FigJ7]the conductance 
distribution P(G) obtained from 1,000,000 configurations 
for a fixed average conductance (G) ~ 3 in the square lat- 
tice [panel (a)] and the honeycomb lattice [panel (b)]. In 
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TABLE II: In square lattice or honeycomb model, corre- 
sponding to eight selected parameter labeled in TabU with 
\Ef\ > A, the average conductance (G) and the second, third, 
ninth moments are listed for the first (with ehD, I and 3 
column) and the second (ehD broken, 2 and 4 column) class 
in the diffusive regime with (G) ~ 3. 




Iar 3.001 1.445 0.929 1.840 1.750 2.196 2.266 2.512 2.645 

2 AR 3.013 1.443 0.923 1.838 1.745 2.184 2.266 2.518 2.660 

4 A r 2.986 1.366 0.862 1.737 1.639 2.061 2.129 2.367 2.486 

6ar 2.976 1.365 0.864 1.736 1.639 2.061 2.129 2.367 2.486 




Iarr 2.998 1.417 0.899 1.805 1.708 2.143 2.218 2.463 2.591 
2 A rr 2.988 1.421 0.901 1.809 1.713 2.148 2.224 2.470 2.602 
5arr 3.005 1.347 0.836 1.713 1.606 2.031 2.091 2.328 2.441 
4 SA r 3.009 1.344 0.833 1.710 1.602 2.027 2.087 2.324 2.438 



the absence of magnetic filed, where 'T' denotes trans- 
pose. Eq.(|10p can be expanded in power series which 
gives multiple Andreev reflections. For qualitative un- 
derstanding, we can focus on the first term in the series, 
i.e., m(e) = t 12 {e)t\ 2 (-e) and T% ] = T 12 (e)T 12 (-e). It 
is similar for the higher order of Ta- Now it is clear 
why we obtain two universal conductance distributions 
for Andreev conductance. For e = (with ehD) the total 
Andreev reflection coefficient T4 is expressed in terms 
of only one type of normal transmission coefficient T(0). 
For e ^ (ehD broken), however, Ta consists of two kinds 
of transmission coefficient T(e) and T(— e) that have the 
completely different statistics. It is the statistical inter- 
ference of T(e) and T(— e) that leads to the new universal 
conductance distribution. 

It should be noted in order to get the uniform statis- 
tical interference, T(e) and T(—e) must be separated far 
enough from each other, i.e., e is larger than Thouless en- 
ergy. The incident energy e (related to condensed energy, 
equal to Ep in our calculation) is so large that it is com- 
parable to energy gap A, so we must go beyond Andreev 
approximation (AA). While in the present works, AA are 
widely used, it is why the present works can't present this 
new symmetry class. We will show [FigJS] in the AA, the 
conductance distribution is smoothly changed with Ep, 
in stead of the two universal functions corresponding to 
Ep = and Ep 7^ in the case with non- Andreev ap- 
proximation (NAA). 



B. Statistical properties in the localized regime 

As we have shown, different universal conductance dis- 
tributions corresponding to Ep — and Ej? / are 
found in the diffusive regime. It has been demonstrated 
numerically-^ that the conductance distribution for a fixed 
(G) in the localized regime seems to be a universal func- 
tion which does not depend on dimensionality (quasi-lD, 
2D and quantum dot systems) and ensemble symmetry 
(COE, CUE or CSE). For normal-superconductor hybrid 
systems, it is interesting to know whether this conclusion 
is still valid. 

There are two ways to examine the universal conduc- 
tance distribution P{G): (1). plot P(G) at each (G) for 
different system parameters to see whether all P(G) col- 
lapse into a single curve. One can only plot P(G) at a few 
selected (G). (2). plot the moments of P(G) as a func- 
tion of (G) to see the universal behavior. However one 
can only plot several moments of conductance. Here we 
focus on the higher order moments 113 and ^,4. In FigJSl 
we plot ^/7i3 [panel(a)] and tfjTi [panel(4)] vs (G) for 2D 
and quasi- ID systems on square and honeycomb lattices. 
Symbols (l)-(9) are described as in panel (b) and labeled 
in TablU From the figure, it is clear that the data do not 
collapse into a single curve. In this calculation, we have 
used only 10,000 configurations per data point which is 
not enough to resolve the universality class if any. To 
improve this, we fix the average conductance (G) and 



this figure, we choose eight parameters from TabU with 
\Ep\ < A. We see that for both square and honeycomb 
lattices, the conductance distributions corresponding to 
Ep = and Ep ^ are clearly different. In addition, for 
each case, Ep — or Ep 7^ 0, conductance distributions 
for square and honeycomb lattices are almost the same, 
as can be seen from TablTTl in which the second, third, 
ninth moments are listed for the parameters labeled 
in TabU corresponding to the first (Ep = with ehD) 
and the second (Ep 7^ where ehD is broken) classes 
with fixed (G) ~ 3. Here, the n-th moment is defined 
as [i n = ([G- (G)]™). In TabHH the n-th moments la- 
beled by "Iar" and u 2ar" correspond to the first class 
in square lattice, they are close to the n-th moments la- 
beled by 'Iarr' and '2,4^' in honeycomb lattice. Since 
the universal behavior is determined only by the sym- 
metry and dimensionality, why there are two universal 
curves for AR? This can be qualitatively understood as 
follows. 

When the energy of incoming electron is within the 
superconducting energy gap, only AR exists. The AR 
amplitude of total NS system Ta is contributed by mul- 
tiple Andreev reflections and can be expressed in terms 
of transmission amplitude t and r in the absence of su- 
perconducting leads and the pure AR matrix ta of the 
only NS interface [not consider the clean or disordered 
normal scattering region] in the following form 3 

T A (e) = Tr[m(e)m+(e)] (9) 

with 

m(e)=t e 12 (e)Mt^(-e) 

M = [I- ^( e )r 2 ^(- e )rl^ T ( e )^ 2 ( e )]^<( e XlO) 

where we have used the electron-hole symmetry relation 
*2x( e ) = *2i*(- e )> r A( e ) = r A h,T i e ) and the symmetry 
relation of normal transmission matrix i 21 (e) = ^i2 T ( e ) m 
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FIG. 8: (Color online) the skewness 71 and the kurtosis 72 vs 
(G) for the ID square or honeycomb lattice with N c = 40 and 
2D square or honeycomb lattice with N c — 40 and N c — 60. 
Different symbols (l)-(9) are described as in panel (b) and 
labeled in TabU] 

TABLE III: Same to TablTTl except we consider localized 
regime with fixed average conductance (G) ~ 0.3. 

Sq (G) ygg ^3 \/M4 jfJZE -g/^e J/Jl7 -y/MS \ffM) 

Ur .3005 0.669 0.986 1.287 1.527 1.728 1.900 2.052 2.192 

2ar .2997 0.669 0.986 1.287 1.527 1.727 1.898 2.049 2.186 

4 AR .3002 0.632 0.936 1.229 1.464 1.658 1.823 1.964 2.087 

6ar .2990 0.631 0.936 1.229 1.464 1.660 1.826 1.970 2.099 

he (G) y/JJq yjMi fyjH y/M6 l/Wi -y/MS %ftM) 

Iarr .2998 0.667 0.982 1.281 1.520 1.718 1.887 2.036 2.170 

2arr -3005 0.669 0.985 1.286 1.526 1.725 1.895 2.044 2.179 

5arr .3001 0.631 0.934 1.226 1.460 1.654 1.817 1.958 2.081 

4sar .3002 0.631 0.936 1.229 1.463 1.658 1.823 1.966 2.094 



calculate higher moments by averaging over 1,000,000 
configurations. In Tab lllll we choose the same set of 
parameters as used in the diffusive regime [TabHTj. and 
tabulate the average conductance (G) and the second, 
third, ninth moments for the fixed (G) ~ 0.3. Similar 
to Tab|nl two universality classes can be identified. The 
first universality class has ehD and consists of data points 
from four different set of parameters labeled by "Iar" 
and "2a/?" (square lattice) and labeled by "Iarr" and 
"2 arr" (honeycomb lattice). The rest of data form the 
second universality class where ehD is broken. Hence it is 
expected that the conductance distributions for Ep = 
(with ehD) and Ep ^ (without ehD) belong to differ- 
ent universality class in the localized regime. This indeed 
can be seen from Figj9ja) and (b) where we have plot- 
ted the conductance distributions of log 10 (G) for Ep = 
and Ep 0. Figl^a) shows the conductance distribution 
with ehD for six different sets of parameters where two of 



-5 -4 -3 -2 -1 r 

l°9 10 (G) 

FIG. 9: (Color online) In the localized regime, corresponding 
to selected parameters labeled in TabUwith \Ef\ = and < 
\Ef\ < A, the conductance distribution P[log 10 (G)] obtained 
from 1,000,000 configurations are plotted in panel (a) and 
panel(b) respectively for the fixed (G) ~ 0.3 in square lattice 
[marked with sq] and honeycomb lattice [marked with he]. In 
addition, we also plot P[log 10 (G)] within AA for Ep = and 
Ep 7^ in panel(a) and panel(b), respectively. 



TABLE IV: Beyond the superconducting Gap A, the average 
conductance (G) and the second, third, ninth moments 
of normal conductance Gjv are listed in the localized regime 
with (G) ~ 0.3. 

NT (G) ^/Jl2 tyjlj ^/JH jfjls ^/JUs y/Wl y/M8 y^9 

lsq .2997 0.311 0.356 0.470 0.546 0.622 0.691 0.755 0.814 

3 sq .2999 0.310 0.355 0.467 0.544 0.620 0.689 0.752 0.812 

l hc -2998 0.310 0.357 0.472 0.551 0.631 0.703 0.771 0.835 

6 hc .3004 0.310 0.351 0.462 0.537 0.612 0.679 0.742 0.801 



them are for AA and the other four are NAA. Obviously, 
they fall into the same universality class. In Fig|9]Jb), we 
show the data for the case with broken ehD. We see that 
four set of data with NAA collapse into a single curve 
indicating the universal conductance distribution that is 
clearly different from Figj9ja). When AA is made, how- 
ever, the conductance distribution depends on Ep which 
is non-universal. The results from FigjS] show that even 
in the localized regime, the Andreev conductance distri- 
butions for Ep — (with ehD) and Ep ^ (ehD broken) 
belong to different universality class. 



C. statistics beyond superconducting gap 

In previous sub-sections, we have studied the statis- 
tical properties of pure AR assisted conductance with 
incident energy \Ep\ < A. In this sub-section, we will 
focus on the case in which the incident energy Ep is 



9 



0.8 

■E 0.6 
"at 

— 0.4 

E 

a 02 

n 

> 0.0 
£" 0.06 

I 0.03 
& 
n 

> 0.00 

0.2 



> 0.0 

0. 



(a1) 



(bi> 



□ 1 B 1 

o 2 6 2 

A3 A3 



■'"/! 



-(d) 



□ 1 

O 2 
A 3 



(a2) 




(b2) 



he 



(c2) 



1 
2 
3 
4 
5 

6 ,il 




01 0.1 1 10 

<G> (e 2 /h) 



0.01 0.1 1 
<G> (e 2 /h) 



10 



FIG. 10: (Color online) Var(G)Norm, Var(G)Andr and 
Var(G) C ross, the three compositions of variance of G vs (G) 
for the 2D square lattice [the left column] and 2D honeycomb 
lattice [the right column] with N c — 40 [open symbols] and 
N c — 60 [symbols with '-']. 




FIG. 11: (Color online) Panel (a): Var(G) Nor m vs (G) for 
the ID [with iV c = 40, corresponding to the crossed symbols] 
or 2D [with 7V C = 40, corresponding to the open symbols and 
N c — 60, corresponding to the symbols with '-'] square lattice 
[symbols (1), (2) and (3)] or honeycomb lattice [symbols (4), 
(5) and (6)]. Panel (b): in a 2D square or honeycomb lattice 
system with N c = 40, corresponding to four selected param- 
eter with Ef beyond A and labeled in TabU conductance 
distribution P[log 10 (G)] exported from 1,000,000 configura- 
tions is plotted. In comparation, we also plot P[log 10 (G)] for 
different parameters in the normal system with A = from 
1,000,000 configurations. 



above A. In this case, conductance is contributed by 
both normal transmission and Andreev reflection. The 
conductance variance Var(G) consists of three terms, the 
Andreev conductance fluctuation Var(G)Andr , the normal 
conductance fluctuation Var(G)Norm and the cross term 
between them Var(G) cross [see Eq.(JBJ) and Eq.flSJ]. In 
FigHUl we plot Var(G) No rm, Var(G) An dr and Var(G) cross 
vs (G) for the 2D square lattice [left panels] and 2D hon- 
eycomb lattice [right panels] with N c = 40 [open sym- 
bols] and N c = 60 [symbols with '-']. Our results can 
be summarized as follows. (1) The Andreev related vari- 
ance Var(G)Andr is drastically suppressed by the disor- 
der. In localized regime [(G) < 1], due to strong dis- 
order, it is completely suppressed to almost zero. As a 
result only Var(G)N orm plays a dominant pole in the lo- 
calized regime. (2) in the localized regime, the dominant 
Var(G)Norm exhibits a universal behavior, i.e., it is in- 
dependent of system parameters (such as Ep, Eq, N c , 
A and so on). In FigHTTa). we plot Var(G)Norm of 2D 
system [Fig|l0Tbl).(b2)] and quasi-lD system for square 
lattice and honeycomb lattice. We find that Var(G)Norm 
in the localized regime is also independent of dimension- 
ality and type of lattice. It is not surprising since in lo- 
calized regime, all AR related process are suppressed by 
the strong disorder. In absence of electron-hole conver- 
sion, statistics of NS system are same as that of normal 
system. 



In order to improve the accuracy in the calculation, 
we also calculate the higher order moments and conduc- 
tance distribution by averaging over 1,000,000 configu- 
rations and tabulate average conductance (G) and the 
second, third, ninth moments for the fixed (G) ~ 0.3 
in Tab II VI It is found that the n-th moment for the 
square lattice and the honeycomb lattice are the same. 
Correspondingly, in Fig lllf b). we plot the conductance 
distribution of log 10 (G) i n a 2D square and honeycomb 
lattices with A = and A ^ 0. The symbols for A ^ 
are labeled as in Tab [I] and the symbols for A = is 
described in FigfTlTb). We see that those data labeled 
by "Ijvt" belong to the first class (Ep = 0), and the 
other data belong to the second class (Ep ^ 0). We can 
see that when the incident energy is above the super- 
conducting gap A, the conductance distributions of NS 
system are almost indistinguishable from that of normal 
system with A = 0. This again confirms that the nor- 
mal transmission is dominant, electron- hole conversion 
and consequently the ehD is irrelevant in the localized 
regime. 



On experimental side, conductance fluctuatio n 3 ^ 33 
and magnetoconductance fluctuation 3 ^ has been mea- 
sured for mono and multi-layer graphene normal systems. 
The conductance fluctuations of normal-superconducting 
hybrid systems (non-graphene) has also been studied^ 
Hence, our results can be checked experimentally. 
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IV. CONCLUSION 

Using the tight-binding model, we have carried out a 
theoretical study on the sample to sample fluctuation in 
transport properties of phase coherent systems with con- 
ventional NS hybrid systems or graphene based NS hy- 
brid systems. Extensive numerical simulations on quasi- 
1D or 2D systems show that (1). When Ep < A, the 
UCF due to AR is found to be roughly doubled compar- 
ing to the system in the absence of the superconducting 
lead. Denoting the increase factor «o through the rela- 
tion rms(Gjvs) = oco rms(Gjv), we found that the differ- 
ence between ao in 2D system and quasi-lD system is 
quite large while the difference is small between Fermi 
electrons and Dirac electrons. (2). Our results show that 
ehD in the NS hybrid system can lead to a new univer- 
sality class. In the diffusive regime we found two slightly 
separated UCF plateaus, one corresponds to the com- 
plete electron- hole symmetry class (with ehD) and the 
other to conventional electron-hole conversion (with ehD 
broken). In addition, the AR conductance distribution 
for the fixed average conductance (G) in diffusive regime 
also confirms that the new universality class can be clas- 



sified using ehD. (3). In the localized regime, we found 
that the conductance distribution is a universal function 
that depends only on the average conductance and the 
ehD. We emphasize that one has to go beyond AA to 
make sure that the AR conductance distribution is uni- 
versal in the localized regime. (4). Finally, when Ep is 
beyond A, normal transport is present. In general, the 
conductance distributions of NS systems and normal sys- 
tems are different. In the localized regime, however, the 
AR is suppressed significantly by the disorder. Hence in 
the localized regime normal transmission dominates the 
transport processes. In this case, the ehD is irrelevant 
and the conductance distribution is a universal function 
that depends only on the average conductance in the lo- 
calized regime. 
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